/**
 * 'MedusaTransformMFCCAdapter'
 * 	
 * - original from the CoMIRVA project
 * - has a dependency on medusa.mobile.client.math.Matrix 
 * and MedusaTransformFFTAdapter
 * 
 * @modified : Nov. 4th, 2011
 * @author : Bin Liu (binliu@usc.edu)
 * 
 */

package medusa.mobile.client;

import java.io.IOException;
import medusa.mobile.client.math.*;

import android.util.Log;

public class MedusaTransformMFCCAdapter {
	
	protected int windowSize;
	protected int hopSize;
	protected float sampleRate;
	protected double baseFreq;

	// fields concerning the mel filter banks
	protected double minFreq;
	protected double maxFreq;
	protected int numberFilters;

	// fields concerning the MFCCs settings
	protected int numberCoefficients;
	protected boolean useFirstCoefficient;

	// implementation details
	@SuppressWarnings("unused")
	private double[] inputData;
	private double[] buffer;
	private Matrix dctMatrix;
	private Matrix melFilterBanks;
	private MedusaTransformFFTAdapter normalizedPowerFFT;
	private double scale;

	/**
	 * Creates a new MFCC object with default window size of 512 for the given
	 * sample rate. The overlap of the windows is fixed at 50 percent. The
	 * number of coefficients is set to 20 and the first coefficient is in use.
	 * The 40 mel-filters are place in the range from 20 to 16000 Hz.
	 * 
	 * @param sampleRate
	 *            float samples per second, must be greater than zero; not
	 *            whole-numbered values get rounded
	 * @throws IllegalArgumentException
	 *             raised if method contract is violated
	 */
	public MedusaTransformMFCCAdapter(float sampleRate)
			throws IllegalArgumentException {
		this(sampleRate, 512, 20, true, 20.0, 16000.0, 40, 256);
	}

	/**
	 * Creates a new MFCC object. 40 mel-filters are place in the range from 20
	 * to 16000 Hz.
	 * 
	 * @param sampleRate
	 *            float samples per second, must be greater than zero; not
	 *            whole-numbered values get rounded
	 * @param windowSize
	 *            int size of window; must be 2^n and at least 32
	 * @param numberCoefficients
	 *            int must be grate or equal to 1 and smaller than the number of
	 *            filters
	 * @param useFirstCoefficient
	 *            boolean indicates whether the first coefficient of the dct
	 *            process should be used in the mfcc feature vector or not
	 * @throws IllegalArgumentException
	 *             raised if method contract is violated
	 */
	public MedusaTransformMFCCAdapter(float sampleRate, int windowSize,
			int numberCoefficients, boolean useFirstCoefficient)
			throws IllegalArgumentException {
		this(sampleRate, windowSize, numberCoefficients, useFirstCoefficient,
				20.0, 16000.0, 40, 256);
	}

	/**
	 * Creates a new MFCC object. 40 mel-filters are place in the range from 20
	 * to 16000 Hz.
	 * 
	 * @param sampleRate
	 *            float samples per second, must be greater than zero; none
	 *            integer values get rounded
	 * @param windowSize
	 *            int size of window; must be 2^n and at least 32
	 * @param numberCoefficients
	 *            int must be grate or equal to 1 and smaller than the number of
	 *            filters
	 * @param useFirstCoefficient
	 *            boolean indicates whether the first coefficient of the dct
	 *            process should be used in the mfcc feature vector or not
	 * @param minFreq
	 *            double start of the interval to place the mel-filters in
	 * @param maxFreq
	 *            double end of the interval to place the mel-filters in
	 * @param numberFilters
	 *            int number of mel-filters to place in the interval
	 * @throws IllegalArgumentException
	 *             raised if method contract is violated
	 */
	public MedusaTransformMFCCAdapter(float sampleRate, int windowSize,
			int numberCoefficients, boolean useFirstCoefficient,
			double minFreq, double maxFreq, int numberFilters, int overlapSize)
			throws IllegalArgumentException {
		// check for correct window size
		if (windowSize < 32) {
			throw new IllegalArgumentException(
					"window size must be at least 32");
		} else {
			int i = 32;
			while (i < windowSize && i < Integer.MAX_VALUE)
				i = i << 1;

			if (i != windowSize)
				throw new IllegalArgumentException("window size must be 2^n");
		}

		// check sample rate
		sampleRate = Math.round(sampleRate);
		if (sampleRate < 1)
			throw new IllegalArgumentException("sample rate must be at least 1");

		// check numberFilters
		if (numberFilters < 2 || numberFilters > (windowSize / 2) + 1)
			throw new IllegalArgumentException(
					"number filters must be at least 2 and smaller than the nyquist frequency");

		// check numberCoefficients
		if (numberCoefficients < 1 || numberCoefficients >= numberFilters)
			throw new IllegalArgumentException(
					"the number of coefficients must be greater or equal to 1 and samller than the number of filters");

		// check minFreq/maxFreq
		if (minFreq <= 0 || minFreq > maxFreq || maxFreq > 88200.0f)
			throw new IllegalArgumentException(
					"the min. frequency must be greater 0 smaller than the max. frequency, which must be smaller than 88200.0");
		;

		this.sampleRate = sampleRate;
		this.windowSize = windowSize;
		this.hopSize = windowSize - overlapSize; // hopSize
		this.baseFreq = sampleRate / windowSize;

		this.numberCoefficients = numberCoefficients;
		this.useFirstCoefficient = useFirstCoefficient;

		this.minFreq = minFreq;
		this.maxFreq = maxFreq;
		this.numberFilters = numberFilters;

		// create buffers
		inputData = new double[windowSize];
		buffer = new double[windowSize];

		// store filter weights and DCT matrix due to performance reason
		melFilterBanks = getMelFilterBanks();
		dctMatrix = getDCTMatrix();

		// create power fft object
		normalizedPowerFFT = new MedusaTransformFFTAdapter(
				MedusaTransformFFTAdapter.FFT_NORMALIZED_POWER, windowSize,
				MedusaTransformFFTAdapter.WND_HANNING);
		// normalizedPowerFFT = new FFT(FFT.FFT_NORMALIZED_POWER, windowSize,
		// FFT.WND_HAMMING);

		// compute rescale factor to rescale and normalize at once (default is
		// 96dB = 2^16)
		scale = (Math.pow(10, 96 / 20));
	}

	/**
	 * Returns the boundaries (start, center, end) of a given number of
	 * triangular mel filters at linear scale. Mel-filters are triangular
	 * filters on the linear scale with an integral (area) of 1. However they
	 * are placed equidistantly on the mel scale, which is non-linear rather
	 * logarithmic. The minimum linear frequency and the maximum linear
	 * frequency define the mel-scaled interval to equidistantly place the
	 * filters. Since mel-filters overlap, an array is used to efficiently store
	 * the boundaries of a filter. For example you can get the boundaries of the
	 * k-th filter by accessing the returned array as follows:
	 * 
	 * leftBoundary = boundaries[k-1]; center = boundaries[k]; rightBoundary =
	 * boundaries[k+1];
	 * 
	 * @param minFreq
	 *            double frequency used for the left boundary of the first
	 *            filter
	 * @param maxFreq
	 *            double frequency used for the right boundary of the last
	 *            filter
	 * @param numberFilters
	 *            int number of filters to place within the interval [minFreq,
	 *            maxFreq]
	 * @return double[] array holding the boundaries
	 */
	private double[] getMelFilterBankBoundaries(double minFreq, double maxFreq,
			int numberFilters) {
		// create return array
		double[] centers = new double[numberFilters + 2];
		double maxFreqMel, minFreqMel, deltaFreqMel, nextCenterMel;

		// compute mel min./max. frequency
		maxFreqMel = linToMelFreq(maxFreq);
		minFreqMel = linToMelFreq(minFreq);
		deltaFreqMel = (maxFreqMel - minFreqMel) / (numberFilters + 1);

		// create (numberFilters + 2) equidistant points for the triangles
		nextCenterMel = minFreqMel;
		for (int i = 0; i < centers.length; i++) {
			// transform the points back to linear scale
			centers[i] = melToLinFreq(nextCenterMel);
			nextCenterMel += deltaFreqMel;
		}

		// ajust boundaries to exactly fit the given min./max. frequency
		centers[0] = minFreq;
		centers[numberFilters + 1] = maxFreq;

		return centers;
	}

	/**
	 * This method creates a matrix containing <code>numberFilters</code>
	 * mel-filters. Each filter is represented by one row of this matrix. Thus
	 * all the filters can be applied at once by a simple matrix multiplication.
	 * 
	 * @return Matrix a matrix containing the filter banks
	 */
	private Matrix getMelFilterBanks() {
		// get boundaries of the different filters
		double[] boundaries = getMelFilterBankBoundaries(minFreq, maxFreq,
				numberFilters);

		// ignore filters outside of spectrum
		for (int i = 1; i < boundaries.length - 1; i++) {
			if (boundaries[i] > sampleRate / 2) {
				numberFilters = i - 1;
				break;
			}
		}

		// create the filter bank matrix
		double[][] matrix = new double[numberFilters][];

		// fill each row of the filter bank matrix with one triangular mel
		// filter
		for (int i = 1; i <= numberFilters; i++) {
			double[] filter = new double[(windowSize / 2) + 1];

			// for each frequency of the fft
			for (int j = 0; j < filter.length; j++) {
				// compute the filter weight of the current triangular mel
				// filter
				double freq = baseFreq * j;
				filter[j] = getMelFilterWeight(i, freq, boundaries);
			}

			// add the computed mel filter to the filter bank
			matrix[i - 1] = filter;
		}

		// return the filter bank
		return new Matrix(matrix, numberFilters, (windowSize / 2) + 1);
	}

	/**
	 * Returns the filter weight of a given mel filter at a given frequency.
	 * Mel-filters are triangular filters on the linear scale with an integral
	 * (area) of 1. However they are placed equidistantly on the mel scale,
	 * which is non-linear rather logarithmic. Consequently there are lots of
	 * high, thin filters at start of the linear scale and rather few and flat
	 * filters at the end of the linear scale. Since the start-, center- and
	 * end-points of the triangular mel-filters on the linear scale are known,
	 * the weights are computed using linear interpolation.
	 * 
	 * @param filterBank
	 *            int the number of the mel-filter, used to extract the
	 *            boundaries of the filter from the array
	 * @param freq
	 *            double the frequency, at which the filter weight should be
	 *            returned
	 * @param boundaries
	 *            double[] an array containing all the boundaries
	 * @return double the filter weight
	 */
	private double getMelFilterWeight(int filterBank, double freq,
			double[] boundaries) {
		// for most frequencies the filter weight is 0
		double result = 0;

		// compute start- , center- and endpoint as well as the height of the
		// filter
		double start = boundaries[filterBank - 1];
		double center = boundaries[filterBank];
		double end = boundaries[filterBank + 1];
		double height = 2.0d / (end - start);

		// is the frequency within the triangular part of the filter
		if (freq >= start && freq <= end) {
			// depending on frequency position within the triangle
			if (freq < center) {
				// ...use a ascending linear function
				result = (freq - start) * (height / (center - start));
			} else {
				// ..use a descending linear function
				result = height
						+ ((freq - center) * (-height / (end - center)));
			}
		}

		return result;
	}

	/**
	 * Compute mel frequency from linear frequency.
	 * 
	 * @param inputFreq
	 *            the input frequency in linear scale
	 * @return the frequency in a mel scale
	 */
	private double linToMelFreq(double inputFreq) {
		return (2595.0 * (Math.log(1.0 + inputFreq / 700.0) / Math.log(10.0)));
	}

	/**
	 * Compute linear frequency from mel frequency.
	 * 
	 * @param inputFreq
	 *            the input frequency in mel scale
	 * @return the frequency in a linear scale
	 */
	private double melToLinFreq(double inputFreq) {
		return (700.0 * (Math.pow(10.0, (inputFreq / 2595.0)) - 1.0));
	}

	/**
	 * Generates the DCT matrix for the known number of filters (input vector)
	 * and for the known number of used coefficients (output vector). Therefore
	 * the DCT matrix has the dimensions (numberCoefficients x numberFilters).
	 * If useFirstCoefficient is set to false the matrix dimensions are
	 * (numberCoefficients-1 x numberFilters). This matrix is a submatrix of the
	 * full matrix. Only the first row is missing.
	 * 
	 * @return Matrix the appropriate DCT matrix
	 */
	private Matrix getDCTMatrix() {
		// compute constants
		double k = Math.PI / numberFilters;
		double w1 = 1.0 / (Math.sqrt(numberFilters));// 1.0/(Math.sqrt(numberFilters/2));
		double w2 = Math.sqrt(2.0 / numberFilters);// Math.sqrt(2.0/numberFilters)*(Math.sqrt(2.0)/2.0);

		// create new matrix
		Matrix matrix = new Matrix(numberCoefficients, numberFilters);

		// generate dct matrix
		for (int i = 0; i < numberCoefficients; i++) {
			for (int j = 0; j < numberFilters; j++) {
				if (i == 0)
					matrix.set(i, j, w1 * Math.cos(k * i * (j + 0.5d)));
				else
					matrix.set(i, j, w2 * Math.cos(k * i * (j + 0.5d)));
			}
		}

		// ajust index if we are using first coefficient
		if (!useFirstCoefficient)
			matrix = matrix.getMatrix(1, numberCoefficients - 1, 0,
					numberFilters - 1);

		return matrix;
	}

	/**
	 * Performs the transformation of the input data to MFCCs. This is done by
	 * splitting the given data into windows and processing each of these
	 * windows with processWindow().
	 * 
	 * @param input
	 *            double[] input data is an array of samples, must be a multiple
	 *            of the hop size, must not be a null value
	 * @return double[][] an array of arrays contains a double array of Sone
	 *         value for each window
	 * @throws IOException
	 *             if there are any problems regarding the inputstream
	 * @throws IllegalArgumentException
	 *             raised if method contract is violated
	 */
	public double[][] process(double[] input) {
		// check for null
		if (input == null) {
			Log.e("MedusaTransformMFCCAdapter", "input data must not be a null value");		
			return null;
		}

		// check for correct array length
		if ((input.length % hopSize) != 0) {
			Log.e("MedusaTransformMFCCAdapter", "Input data must be multiple of hop size (windowSize/2).");
			return null;
		}

		// create return array with appropriate size
		double[][] mfcc = null;

		if (windowSize > hopSize) {
			mfcc = new double[(input.length / hopSize) - 1][numberCoefficients];
		} else if (windowSize == hopSize) {
			mfcc = new double[input.length / hopSize][numberCoefficients];
		}

		// process each window of this audio segment
		for (int i = 0, pos = 0; pos < input.length - hopSize; i++, pos += hopSize)
			mfcc[i] = processWindow(input, pos);

		if (windowSize == hopSize) {
			mfcc[mfcc.length - 1] = processWindow(input, input.length - hopSize);
		}

		return mfcc;
	}

	/**
	 * Returns the window size.
	 * 
	 * @return int the window size in samples
	 */
	public int getWindowSize() {
		return windowSize;
	}

	/**
	 * Transforms one window of MFCCs. The following steps are performed: <br>
	 * <br>
	 * (1) normalized power fft with hanning window function<br>
	 * (2) convert to Mel scale by applying a mel filter bank<br>
	 * (3) Conversion to db<br>
	 * (4) finally a DCT is performed to get the mfcc<br>
	 * <br>
	 * This process is mathematical identical with the process described in [1].
	 * 
	 * @param window
	 *            double[] data to be converted, must contain enough data for
	 *            one window
	 * @param start
	 *            int start index of the window data
	 * @return double[] the window representation in Sone
	 * @throws IllegalArgumentException
	 *             raised if method contract is violated
	 */
	public double[] processWindow(double[] window, int start)
			throws IllegalArgumentException {
		// number of unique coefficients, and the rest are symmetrically
		// redundant
		int fftSize = (windowSize / 2) + 1;

		// check start
		if (start < 0)
			throw new IllegalArgumentException("start must be a positve value");

		// check window size
		if (window == null || window.length - start < windowSize)
			throw new IllegalArgumentException(
					"the given data array must not be a null value and must contain data for one window");

		Boolean cont = false;
		// just copy to buffer and rescaled the input samples according to the
		// original matlab implementation to 96dB
		for (int j = 0; j < windowSize; j++) {
			buffer[j] = window[j + start];

			if (buffer[j] != 0) {
				buffer[j] *= scale;
				if (!cont)
					cont = true;
			}

		}

		if (cont) {
			// perform power fft
			normalizedPowerFFT.transform(buffer, null);

			// use all coefficient up to the nyquist frequency
			// (ceil((fftSize+1)/2))
			Matrix x = new Matrix(buffer, windowSize);
			x = x.getMatrix(0, fftSize - 1, 0, 0); // fftSize-1 is the index of
													// the nyquist frequency

			// apply mel filter banks
			x = melFilterBanks.times(x);

			// to db
			double log10 = 10 * (1 / Math.log(10)); // log for base 10 and scale
													// by factor 10
			x.thrunkAtLowerBoundary(1);
			x.logEquals();
			x.timesEquals(log10);

			// compute DCT
			x = dctMatrix.times(x);

			return x.getColumnPackedCopy();
		} else {
			return new double[numberCoefficients];
		}
	}

}
